Genetic Characteristics and Phylogeographic Dynamics of Lagoviruses, 1988–2021

Rabbit haemorrhagic disease virus (RHDV), European brown hare syndrome virus (EBHSV), rabbit calicivirus (RCV), and hare calicivirus (HaCV) belong to the genus Lagovirus of the Caliciviridae family that causes severe diseases in rabbits and several hare (Lepus) species. Previously, Lagoviruses were classified into two genogroups, e.g., GI (RHDVs and RCVs) and GII (EBHSV and HaCV) based on partial genomes, e.g., VP60 coding sequences. Herein, we provide a robust phylogenetic classification of all the Lagovirus strains based on full-length genomes, grouping all the available 240 strains identified between 1988 and 2021 into four distinct clades, e.g., GI.1 (classical RHDV), GI.2 (RHDV2), HaCV/EBHSV, and RCV, where the GI.1 clade is further classified into four (GI.1a–d) and GI.2 into six sub-clades (GI.2a–f). Moreover, the phylogeographic analysis revealed that the EBHSV and HaCV strains share their ancestor with the GI.1, while the RCV shares with the GI.2. In addition, all 2020–2021 RHDV2 outbreak strains in the USA are connected to the strains from Canada and Germany, while RHDV strains isolated in Australia are connected with the USA-Germany haplotype RHDV strain. Furthermore, we identified six recombination events in the VP60, VP10, and RNA-dependent RNA polymerase (RdRp) coding regions using the full-length genomes. The amino acid variability analysis showed that the variability index exceeded the threshold of 1.00 in the ORF1-encoded polyprotein and ORF2-encoded VP10 protein, respectively, indicating significant amino acid drift with the emergence of new strains. The current study is an update of the phylogenetic and phylogeographic information of Lagoviruses that may be used to map the evolutionary history and provide hints for the genetic basis of their emergence and re-emergence.


Introduction
Lagoviruses make up a separate genus, e.g., Lagovirus, in the family Caliciviridae and are the causative agents of severe diseases in the European rabbit (Oryctolagus cuniculus) as well as several hare (Lepus) species. During the 1980s, lagoviruses were first reported when a new epidemic disease in the European brown hare population of Sweden began to cause deaths, and the disease was termed European brown hare syndrome (EBHS) [1]. In the following years (during 1984), a similar disease was reported in China among the farmed rabbits [2] and subsequently was identified across the globe [3], devastating the environment and causing significant economic losses [4][5][6][7]. The disease was termed rabbit 1 and 3 days, and rabbits usually succumb within 12 to 36 h after the onset of fever (>40 • C). Clinically, RHD causes hepatic necrosis, respiratory manifestations (tracheitis, dyspnoea, and cyanosis), neurological disorders, and sudden death [8]. Both the EBHSV and RHDV share typical features with other members of the family Caliciviridae. RHDV mature virions are nonenveloped icosahedral virus particles of 30-40 nm in diameter [8], containing a positive-sense, single-stranded RNA genome of~7.4 kb in length, with two slightly overlapping open reading frames (ORFs), ORF1 and ORF2. ORF1 encodes a polyprotein that, during virus replication, splits into the major capsid protein VP60 and seven non-structural proteins, including p16, p23, helicase, p29, VPg, protease, and RNAdependent RNA polymerase (RdRp), while ORF2 encodes the minor structural protein VP10 [28]. A sub-genomic RNA (2.2 kb) is collinear to the 3 end of the virus genome. VP60 is the immunogenic protein and the major structural capsid protein of RHDV [16]. VP60 comprises three domains, the N-terminal arm (NTA), the shell (S), and the protrusion (P), where a short hinge connects the S and P domains [29]. The P domain splits into two sub-domains: P1 and P2. The P2 sub-domain plays an important role in recognizing the histo-blood group antigens (HBGAs) that facilitate the binding of the virus to the host tissue, acting as a receptor [30]. The VP60 gene is widely used to deduce the phylogenetic relationships among Lagovirus isolates [24,28,31].
Phylogenetically, Lagoviruses were classified into the two genogroups GI and GII, where GI consists of RHDVs and RCV, while GII consists of EBHSV and HaCV based on the complete VP60 coding nucleotide sequences [10]. Similarly, a recently proposed taxonomy classified RHDVs into two main genotypes: the classical RHDV named GI.1, which has four further variants: GI.1a to GI.1d [10,32], and the new variant RHDV2 named GI.2, first detected in 2010 in France, causing massive declines in the European rabbit populations [17]. Based on VP60 coding nucleotide sequences, RHDV2 formed a new genetic group, more closely related to classical RHDV and RHDV-related viruses [17]. RHDV2 currently dominates many countries, causes death in young rabbits (<two months old), and infects rabbits of different ages. RHDV2 has a much broader host range, including hares (Lepus spp.) and cottontail rabbits (Sylvilagus spp.) [7]. The large divergence between GI.1 and GI.2 (>15%) might be responsible for the incomplete/low protection of vaccines developed against GI.1 and GI.2 outbreaks [17].
While the main route of RHDV infection is oral transmission [33], the virus is environmentally stable, highly infectious, and transmissible by close contact or by contact with fomites such as contaminated fur, clothing, or cages [33,34]. Indirect arthropod vectors, including blow flies or flesh flies, have also been implicated in the spread of RHDV [34]. The virus spreads worldwide and is considered endemic in most countries, leading to enormous economic loss in the global rabbit industry and impacting human society and natural ecosystems. It has been previously demonstrated that there are many strains of RHDV circulating in different rabbit populations, most of which show distinct epidemiological, pathogenetic, and genetic characteristics [35][36][37][38][39]. Thus, it is important to gain more insight into the genomic evolutionary characteristics of Lagoviruses to help formulate the RHD and EBHS management, virus detection, and exploration of more effective vaccines and treatments. Previously, phylogenetic analyses were carried out based on partial or complete sequences of the VP60 gene. However, with the emergence of genetic variations and recombination between Lagovirus strains, it may not be sufficient to explain the genetic diversity of these viruses and in-depth molecular epidemiology. In this study, we analyzed the phylogenetic relatedness among the RHDV, RCV, EBHSV, and HaCV strains isolated worldwide from 1988 to 2021 using the full-length, complete VP60 and complete VP10 genomic sequences and mapped the phylogeographic network of the viruses. The results might better define the previous evolutionary characteristics and phylogeographic distribution of RHDV, EBHSV, RCV, and HaCV strains and provide helpful hints for the genetic basis of their emergence and re-emergence.

Phylogenetic Tree Construction and Genomic Similarity Analysis
The nucleotide sequences were aligned using ClustalW in the BioEdit version 7.2.5 package [40]. The unrooted Maximum Likelihood (ML) phylogenetic trees were constructed using the IQ-TREE multicore version 1.6.12 with the best-fitting model SYM + I + G4 for the full-length, GTR + F +I + G4 for complete VP60 coding sequences and TIM + F + I + G4 model for the complete VP10 coding sequences, with 1000 bootstraps [41]. Bootstrap values are shown at each node in the phylogenetic tree. FigTree v1.4 was used to visualize and modify the trees (http://tree.bio.ed.ac.uk/software/figtree/, accessed on 22 September 2022). SimPlot v.3.5.1 was used for the genetic similarity map of the sequences [42].

Phylogeographic Network Analysis
The phylogeographic network depicts genetic linkages between the intra-specific sequences and infers relationships for interpreting population genetic data [43]. Thus, all the sequences were analysed for SNPs (single nucleotide polymorphisms) and Haplotypes using the DnaSP v6 [44]. The MSN (Minimum Spinning Network), implemented by PopArt v1.7 [43], was inferred to visualize the phylogeographic network of the strains. The geographic distribution of the 240 full-length sequences isolated during 1988-2021 in nineteen countries was employed for building the network.

Recombination Analysis
Recombination events among 240 full-length nucleotide sequences of Lagoviruses were analyzed using Recombination Detection Program 4 software (RDP4) [45]. Each of the seven algorithms embedded in the RDP4 package (GENECONV, PhylPro, MaxChi, RDP, Bootscan, SiScan, and Chimaera) were used to identify potential recombination events. In this report, a recombination event confirmed by at least four of seven methods was accepted as real.

Amino Acids Variability Analysis
All ORF1 and ORF2 complete nucleotide sequences of 240 Lagovirus strains were retrieved separately from the NCBI database and aligned with ClustalW using the MEGA 11 [46,47]. Nucleotide sequences were translated and edited using BioEdit version 7.2.5 [40]. The amino acid variability was determined using the Protein Variability Server (PVS) with the Wu-Kabat variability coefficient method [48]. The variability coefficient was calculated with the formula variability = N × k/n, where N represents the number of sequences in the alignment, k represents the number of different amino acids at a given position, and n corresponds to the time that the most commonly recognized amino acid at that position is available [49].

Genotyping of Lagoviruses Based on the Full-Length Genomes or Coding Sequences of VP60 and VP10
We assessed a total of 240 full-length genome sequences of Lagoviruses (isolated between 1988-2021) available on the NCBI GenBank database to determine the global variation and phylogenetic and phylogeographic characteristics of RHDV, EBHSV, RCV, and HaCV. Following multiple sequence alignments, the full-length genome-based phylogenetic tree was constructed with the best-fit model SYM + I + G4 based on 1000 bootstrap replications as suggested by IQ-TREE multicore version 1.6.12 [41]. There were a total of 7359 positions in the final dataset of the Maximum Likelihood (ML) method.
Since the Lagovirus genomes showed great diversity, we further analyzed their genetic similarities using fifteen representative full-length strains from each clade and sub-clade and GS/YZ(China-2011-Deer) (GenBank ID: MN478485.1) as the query strain ( Figure 2A,B). In consistency with the phylogenetic tree results, the genetic similarity map showed a wide genetic diversity across the complete Lagovirus genomes (<95%). The results indicated that the RCV, HaCV/EBHSV, and RHDV genomes were separated into different groups, which is consistent with the phylogenetic results ( Figure 1). Furthermore, the RHDV varied into two distinct groups (GI.1 and GI.2) from 1 to 5250 nucleotide positions, encoding p16, p23, RNA-helicase, p29, VPg, peptidase, and RdRp, respectively. Nucleotide positions 5250-7437, encoding VP60 and VP10, revealed the greatest genetic variation indicated by the lowest similarity levels (<70% and <85%, respectively). The p23 and RdRp regions were found to be the most conserved regions of RHDV. However, RCV and HaCV/EBHSV exhibited great diversity across full-length genomes with the lowest similarity (<75%) in the genomic regions nt 1-800 and nt 6000-7000, encoding p16, p23, and VP60, respectively (Figure 2A,B). Therefore, our genomic similarity results speculate the occurrence of recombination between the Lagovirus strains belonging to the four clades (RCV, HaCV/EBHSV, GI.1, and GI.2).
HaCV/EBHSV exhibited great diversity across full-length genomes with the lowest similarity (<75%) in the genomic regions nt 1-800 and nt 6000-7000, encoding p16, p23, and VP60, respectively (Figure 2A,B). Therefore, our genomic similarity results speculate the occurrence of recombination between the Lagovirus strains belonging to the four clades (RCV, HaCV/EBHSV, GI.1, and GI.2). As the VP60 coding region is widely used for classifying RHDVs, we constructed additional phylogenetic trees based on VP60 and VP10 complete coding sequences. Similar to the full-length genome-based genotyping tree results, VP60 phylogenetic analysis identified four distinct clades (RCV, HaCV/EBHSV, GI.1: RHDV, and GI.2: RHDV2) with four sub-clades within the GI.1 clade (GI.1a-d), similar to full-length tree, with four subclades within GI.2 clade (GI.2a-d) (Figure 3, Supplementary Figure S2), unlike the six subclades within the full-length sequence tree (Figure 1, Supplementary Figure S1). In addition, there were significant differences between the full-length genome and VP60 coding sequence-based grouping of the strains. For example, the virus strain GI.1d/O To further explore the phylogenetic relatedness, we also used the available complete VP10 nucleotide sequences (a total of 240) to determine the genotyping of the Lagoviruses (Supplementary Figure S3). The VP10 nucleotide sequence-based tree grouped all the strains into three major clades (GI.1:RHDV, RCV, and GI.2: RHDV2/HaCV/EBHSV). However, the sub-clade genotyping was unclear, especially for the GI.2 clade. Many strains were shown to be grouped into different clades and sub-clades. For example, the virus GI.1d/O cun/FR/2009/09-03 (GenBank ID: MT628290.1, France-2009) grouped into GI.2 (RHDV2) and GI.1b by full-length genome and VP60, respectively, was grouped into GI.1 by VP10; virus strain P95 (GenBank ID: KJ943791.1, Portugal-1996) was grouped into GI.2f (RHDV2) and GI.1d by full-length sequence and VP60, respectively and the strain SD (GenBank ID: Z29514.1, France-1993) grouped into GI.1d by full-length sequence and VP60, respectively, were grouped into possible GI.1c by VP10 (Supplementary Figure S3). Thus, the VP60-based phylogeny is closer to the RHDV full-length genome sequencebased genotyping in contrast to the VP10-based phylogeny, in which the short VP10 coding sequence may provide limited genetic information. Nevertheless, the full-length genomic sequence-based phylogeny provides a more reliable and robust classification of all the Lagoviruses.  -1996) were grouped into GI.2 (GI.2e and GI.2f, respectively) in the full-length genome-based phylogenetic tree, while they were grouped into GI.1b and GI.1d, respectively, in the VP60-based phylogenetic tree. Similarly, the strain GI.2/O cun/FR/2013/13-165 (GenBank ID: MN737112.1, France-2013) was classified into GI.2a (RHDV2) by VP60, while it was put in GI.1d by the full-length genome sequence analysis.

Phylogeographic Analysis of Full-Length Lagovirus Genomes
To further explore the regional level spread of Lagoviruses, we constructed a phylogeographic network of all the full-length genome sequences. The network analysis indicated a great diversity of RHDVs, showing two main network clusters shown as G1.1(RHDV) and G1.2(RHDV2), respectively, with multiple mutational sub-branches (Figure 4). The analysis identified that the FRG-USA strain (GenBank ID: NC_001543.1) genome isolated in 2000 is haplotype of the FRG-Germany strain (GenBank ID: M67473.1) isolated in 1991. Interestingly, all the strains isolated in Australia before 2010, three from Poland, two from New Zealand, and one from Italy, were found to be connected to the FRG-USA/Germany haplotype through short mutational branches. In addition, the strains isolated in Germany were shown to be the most diverse, spreading across most sub-branches, which is consistent with our ML phylogenetic tree (Figure 1, Supplementary Figure S1).   Furthermore, the phylogeographic analysis revealed that all the HaCV/EBHSV strains clustered together and were sharing their ancestor with the G1.1 (RHDV) strains, connected to the RHDV_GER-NW_D51-1.L00911_2014_Germany (GenBank ID: LR899189.1) strain through a short mutational branch. On the other hand, all the RCV strains clustered with the GI.2 (RHDV2) strains that suggests their common ancestor with the G1.2 (RHDV2) strains. The RCV strains are connected through a short mutational branch with the CBMad17-1_Portugal-2017-Rabbit (GenBank ID: MF407655.1) strain, while connected to the RHDV_GER-NW_D61-2.L00910_2014_Germany (GenBank ID: LR899186.1) strain through a long mutational branch. It is worth mentioning that, within the G1.2 (RHDV2) Network-cluster, all 2020-2021 USA RHDV2 outbreak strains formed a sub-branch and were connected to the BC_Canada_WIN-AH-2018-OTH-0029 strain isolated in Canada in 2018 (GenBank ID: MT900571.1); together, they were connected through a long mutational branch to the same RHDV_GER-NW_D61-2.L00910_2014_Germany strain (GenBank ID: LR899186.1) (Figure 4). These results indicate potential genetic exchange among the Lagoviruses that may have occurred by animal trade across the globe.

Recombination Pattern of Lagovirus Full-Length Genome
Since the genomic similarity and phylogeographic analyses indicated possible genetic exchanges and evidence of recombination have been reported for part numbers of Lagoviruses [50][51][52], we used the largest dataset (240 full-length sequences) and the improved methodology to systemically assess the occurrence of recombination among the Lagovirus strains. Recombination patterns and the genomic breakpoints mapping were evaluated using the seven algorithms embedded in the RDP4 software package [45]. We identified a total of six recombination events, two of which were inter-genotype (Events 2 and 3), while the four others were intra-genotype (Events 1 and 4-6) ( Table 2). Among them, five recombination events (Events 1-5) occurred within the RHDV strains, one within the RCV, while no recombination event was detected within or among the HaCV/EBHSV strains. Importantly, four recombination events occurred within the VP60/VP10 coding region, with beginning breakpoints at (5318, 5307, 5341, and 5304 nt) and ending at (7366, 7348, 7366 and 7332 nt), respectively, corroborating the genomic similarity analysis findings ( Figure 2). One Event (F77-3, Event 5) occurred within the major capsid proteinVP60, mapped at (nt 6597, nt 7142), and one (RHDV/GER-NW/EI70-1.L03601/2016, Event 4) within the RdRp region at (nt 3715, nt 5306) ( Figure 5).  HaCV/EBHSV strains. Importantly, four recombination events occurred within the VP60/VP10 coding region, with beginning breakpoints at (5318, 5307, 5341, and 5304 nt) and ending at (7366, 7348, 7366 and 7332 nt), respectively, corroborating the genomic similarity analysis findings (Figure 2). One Event (F77-3, Event 5) occurred within the major capsid proteinVP60, mapped at (nt 6597, nt 7142), and one (RHDV/GER-NW/EI70-1.L03601/2016, Event 4) within the RdRp region at (nt 3715, nt 5306) ( Figure 5).  To provide further positive evidence of the recombination events, we constructed three phylogenetic trees based on three fragments of the virus genome. The first fragment (nt 1-3500) encoded for p16, p23, RNA-helicase, p29, VPg, and peptidase; the second fragment (nt 3501-5200) encoded RdRp; while the third fragment (5201-7437) encoded VP60 and VP10. Recombinant, major, and minor parental sequences used to test the recombination are indicated in the short fragment-based trees in red, orange, and blue, respectively. As shown in Supplementary Figure S4, recombinants in Events 1-3 are nested with their major parents in both nt 1-3500 (Supplementary Figure S4A) and nt 3501-5200 (Supplementary Figure S4B)-based phylogenetic trees, but nested with their minor parents in the nt 5201-7437-based phylogenetic tree (Supplementary Figure S4C). Recombinant in Event 4 is nested with its major parent in the nt 1-3500-based phylogenetic tree (Supplementary Figure S4A) but is nested with its minor parents in the nt 3501-5200based tree (Supplementary Figure S4B). These results suggest that the recombination events in RHDV and RCV may drive the rise of new sub-clades.

Amino Acids Variability Landscape of Lagovirus Proteins
We assessed the amino acid variability landscape across the ORF1-encoded polyprotein and ORF2-encoded VP10 of 240 Lagoviruses using the Wu-Kabat variability coefficient implemented by PVS (protein variability server). The polyprotein consensus sequence contained a total of 2385 amino acids that split into VP60 and seven non-structural proteins, including p16, p23, helicase, p29, VPg, protease, and RdRp during the virus replication. Similarly, the VP10 consensus sequence was 118 amino acids long. The Wu-Kabat variability coefficient showed great variability across the polyprotein, with multiple regions having values greater than the estimation limit (1.00). The highest variable region identified was aa 1-270 (highest recorded value = 47) in the p16 and p23 region ( Figure 6A), followed by aa 670-800 in p29 ( Figure 6A,B) and aa 2101-2251 in VP60 ( Figure 6C). In addition, multiple short regions also indicated high variability, e.g., aa 1153-1157 (highest value = 17), aa 1352-1358 (highest value = 19), aa 1685-1737, and aa 1755-1800. On the other hand, the most variable region of VP10 is from aa 58 to 72 (highest value = 12) ( Figure 6D). These results indicate that the Lagovirus amino acids varied greatly during 1988-2021. On the other hand, the most variable region of VP10 is from aa 58 to 72 (highest value = 12) ( Figure 6D). These results indicate that the Lagovirus amino acids varied greatly during 1988-2021.
Recently, the analysis based on VP60 coding sequences classified the lagoviruses into two major genotypes GI and GII, where all the RHDV strains (RHDV and RHDV2) together with the RCVs were placed within GI, while the EBHSV and HaCV were within the GII [10]. However, the availability and accumulation of new Lagovirus strain sequences require revising the current genotyping for a more explicit phylogenetic system placing the existing and new strains accordingly. To this aim, we evaluated 240 full-length sequences of Lagoviruses available in the NCBI GenBank Database, collected between 1988 and 2021, including 183 RHDV, 46 RCV, 8 EBHSV, and 3 HaCV complete genomes. Our full-length genome sequence-based ML phylogenetic tree showed that all the Lagovirus strains were classified into four major clades, e.g., GI.1 (consisting of classical RHDV), GI.2 (consisting of RHDV2), RCV, and HaCV/EBHSV), with a further four RHDV sub-clades within GI.1 (GI.1a-d) and six within GI.2 (GI.2 a-f). Similarly, the complete VP60-based phylogenetic tree exhibited two clades (GI.1 and GI.2) and four sub-clades within GI.1 (Figure 3), supporting the previous studies' classification for RHDVs [10,32]. The RHDV strains that emerged first in France in 2010 [17] and subsequently spread across Europe have been reported to form a unique branch independent of the classical RHDV that is identified as GI.2 (RHDV2). In our report, RHDV2 (GI.2) is recognized as a distinct clade and can be further divided into six further sub-clades (GI.2a-f) ( Figure 1). Therefore, our analysis results add insights into the evolutionary history of the virus, indicating that the newly emerging strains are forming independent phylogenetic branches in clade GI.  [54]. In contrast, the herein phylogenetic tree displayed that all the America RHDV strains collected between 2020 and 2021 fall into GI.2 clade within the GI.2b sub-clade; furthermore, the North of the USA (NY) and South of the USA (NM, TX, and AZ) strains collected in 2020-2021 are shown clustering within the same sub-genogroup (GI.2b) (Figure 1).
As evidenced by the phylogenetic classification, in the network analysis, RHDV fulllength sequences seem to follow two major network clusters, corresponding to GI.1 and GI.2, each of which indicated multiple sub-branches. Within the first network cluster (GI.1), the FRG-Germany strain isolated in 1991 and the FRG-USA strain isolated in 2000 were shown to have evolved from a common ancestor (haplotype), to which the before 2010 Australia sub-branch is linked with other strains from Poland, New Zealand, and Italy; concomitantly, the classical USA strains (n = 4) are shown appertaining to another sub-branch within the same network cluster. The HaCV/ EBHSV strains are connected to this first cluster of RHDVs (GI.1).
In the second network cluster (GI.2), RHDV strains isolated in France and Germany were found in all branches and sub-branches, probably causing the ongoing emergence of the GI.2 virus. All the RCV strains were linked to the GI.2 cluster. In this network cluster, all RHDV strains of the new 2020-2021 USA outbreak clustered together in a sub-branch genetically close to Canada and Poland strains; all were connected to the Germany RHDV strain (2014), separately from the newly emerging Africa strain (GenBank ID: MW118115.1), that clustered into a different sub-branch and was directly connected to the Netherland NL2016 strain through a short mutational branch (77 mutations). The latter supports our phylogenetic classification and corroborates the Ambagala et al. report demonstrating that the Ghana strain (GenBank ID: MW11811) showed a great resemblance of~98.84% at the nucleotide level to the Netherland NL2016 strain, suggesting the Netherland 2015-2017 RHDV outbreaks as the causative agent of RHDV2 in Ghana [4]. As the origin of GI.2 (RHDV2), RCV, and HaCV/EBHSV is still uncertain, these findings are of great importance to understand the dynamics of viral evolution and dispersal and to promote the control of rabbit populations.
GI.2 was reported for the first time in 2010 as a novel RHDV variant [17] that led to death among vaccinated rabbits [56], revealing more divergence compared to the known Lagoviruses, with increased virulence and pathogenicity noticed after 2016 [57]. Our results identified that strain RHDV-SD (GenBank ID: Z29514.1) collected in France 1993 is a recombinant type (GI.2a/GI.1b) that was circulating even before the first announcement of GI.2 outbreaks. RHDV-SD resulted from recombination between European strains from Portugal (1994) and Germany (2007), respectively.
Such as the other genera of Caliciviridae (Vesivirus, Sapovirus, and Norovirus (NoV)), showing that recombination events occur at the starting points encoding the major VP1 capsid protein [58], recombination in RHDVs has also been suggested [59,60] and was reported within the RHDVb strains (RHDV2), with a single major breakpoint located in the 5 region of VP60 gene [61]. Here, four identified recombination events happened within the major capsid protein VP60 (Events 1-3 and 5), and three involved in the minor capsid protein VP10 (Events 1-3). The identified recombinant Events 2 and 3 seemed to involve a combination of non-structural proteins from GI.1 with major and minor capsid proteins (VP60/VP10) from GI.2, which are of great importance for virus genotyping [28]. VP60 is the major structural and immunogenic protein [62] that is widely used to develop vaccines against rabbit haemorrhagic disease (RHD) by producing viral-like particles (VLPs) with a baculovirus expression vector system (BEVS) [63][64][65]. We identified multiple amino acid variations in the VP60 protein, which should be considered while designing any future vaccines. Thus, our results may provide valued data for evaluating and developing vaccines specified to the VP60 protein. Further, the amino acid variations in the RHDV proteins using the Wu-Kabat coefficient showed great variability throughout the protein sequences (>1.00 threshold), suggesting amino acid genetic drift with the emergence of new RHDV and RCV variants. Therefore, we speculate that substantial genetic exchanges and recombination in the RHDV and RCV viral genomes were involved in generating new RHDV lineages.
Currently, the finding of new species of Lagoviruses, the co-infection of multiple species, and the genomic recombination between different species constitute a complex picture of virus spreading and emerging, which is a challenging threat to the rabbit population and infection control. This study provides an updated phylogenetic and phylogeographic landscape of Lagoviruses based on full-length and individual genomic fragments, which illuminate the global distribution and classification of the circulating strains. In addition, we recommend carefully choosing the RHDV isolates during vaccine development to promote the prevention and control management of RHDVs.

Institutional Review Board Statement:
The authors confirm that the ethical policies of the journal, as noted on the journal's author guidelines page, have been adhered to. No ethical approval was required as the data analyzed in this article have been collected.

Informed Consent Statement: Not applicable.
Data Availability Statement: The data are available on the National Center for Biotechnology Information (NCBI) GenBank Database website (https://www.ncbi.nlm.nih.gov/nucleotide/).